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ABSTRACT 


In this thesis we study several computer implementa- 
tions of the thinning algorithm, a new method for generating 
non-homogeneous Poisson processes. The method, developed 
by Professor P.A.W. Lewis, Naval Postgraduate School, 
Monterey, California, and G.S. Shedler, IBM Research 
mieotacOlry, San Jose, California, is valid for Poisson 
processes with any given intensity function. The basic 
thinning algorithm is modified to exploit several refine- 
ments which reduce computer execution time by approximately 
one-third. The basic and modified thinning programs are 
compared with a previous algorithm of Lewis and Shedler, 
the Poisson decomposition and gap-Statistics algorithm, 
which is easily implemented for Poisson processes with 
Mieensity functions of the form exp(a ta,tta,t*). The 
thinning programs are competitive in both execution time 
and computer memory requirements. One program implementa- 
tion generates the events in a Poisson process one at a 


time; another program implements the algorithmic refinements 


which improve efficiency. 
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Pees ODUCTION 


The Poisson process iS a widely known and studied 
stochastic process. It is frequently used to describe 
random arrivals at some type of service facility such as 
a service station fuel pump or a bank teller's window. 
memes MOSt common form, the "rate" of these arrivals is 
considered to be constant over time. This is the homogeneous 
Poisson process which has the familiar property that times 
between arrivals (or events) are exponentially distributed 
With mean equal to the inverse of the rate. 

The assumption of a constant rate, or homogeneity, is 
at best tenuous when applied to real world data. For exam- 
Mey che rate of arrivals at a traffic light typically 
Meeres from very hich during rush hours to very low in the 
Seeery Morning. In addition to this cyclic time-of-day 
effect, arrival rates may exhibit longer term increases or 
decreases. Further, these effects may be superimposed upon 
shorter term effects to produce a more complex rate which 
Varies with time. These processes for which arrival rates 
Vary with time may often be represented by a non-homogeneous 
Poisson process, that is, a Poisson process with a time 
dependent rate of arrival. 

The generic term of Poisson process includes then both 
homogeneous and non-homogeneous Poisson processes. LEWIS 


and syrpuEeR [Ref. 1] define the Poisson process generally in 


aik@) 





terms of a monotone non-decreaSing right-continous function 


Meeewhich is bounded in any finite interval. Then the 


memmeer Of points, N{t'',t'), in any finite interval has a 
feeeocon G1Stribution with parameter u(t'',t’) = A(t') -— A(t'"'). 
MieyesOr example in (0,t'}], with t’' > 0, P{N(t'',t") =n} 

: n'O0,. 

= PiN, ,=n} ae Vie where ae MeCO ee) Ge Or), 


The right derivative A(t) of A(t) will be assumed to 
feeeeeand 15 called the rate function or intensity function 
of the process. A(t) 1s called the integrated rate function 
Meenas the interpretation that for t > 0, A(t) - A(0) = E[N,]. 
For the homogeneous Poisson process iX(t) iS a constant, 
emoeeA, and thus the integrated rate function is simply the 


mmecuet: OL A and t, 1.e. the expected value of N, = N(0O,t). 


Ec 
While simulation of homogeneous Poisson processes is 
relatively straightforward, the non-homogeneous Poisson 
process iS more problematical. Times between events are 
not exponential in the general case and simulation has 
typically been tailored to specific classes of intensity 
Memmerli1OnsS. LEWIS and sHEpLER [Ref. 1] list three general 
methods for simulating non-homogeneous Poisson processes 
and one method for a special rate function. The general 
methods include the time scale transformation method and 
the conditioning and order-statistics method. The special 
methed is the gap-statistics method, a method which is 


particular to the degree-one exponential polynomial intensity 


mijention, i.e. those of the form A(t) = exp (bp 2° bjt). 


dal 





Implementation of the general methods on a computer may 
pose special problems. Often the inverse of the integrated 
rate function is not explicit and must be computed numerically. 
Other problems in implementation generally result in lower 
efficiency, as measured by execution time or computer storage 
requirements or both. 

One class of intensity functions which is of general 
interest is the degree-two exponential polynomial family. 
meetieis, those with intensity function of the form 


A(t) = exp (ap 15 a,c +a on This family of functions has 


wa 
the property of being positive for all values of t, a 
necessary condition for an intensity function. Additionally, 
by varying the magnitude and sign of the coefficients, the 
exponential polynomial of degree two can be made to be mono- 
tone increasing or decreaSing over time, as well as increasing 
and then decreasing, or vice versa. Use of this type of 
intensity function also leads to Statistical procedures 

which are relatively simple. 

LEWIS and SHEDLER [Ref. 2] proposed a new method of 
generating the non-homogeneous Poisson process with degree- 
two exponential polynomial intensity function. It involves 
decomposition of the degree-two exponential polynomial 
Mm@eenmsity function, A(t), into two functions, a degree-one 


exponential polynomial function, (t), and a difference 


Ay, 
maaction, Ay (t) = A(t) - A, (t). This procedure allows the 


points in the degree-one exponential polynomial event stream 


eZ 





to be generated using the gap-statistics method, which is 
highly efficient when implemented on a computer. The 
remaining points with intensity function hy ft) are generated 
by other methods and then merged with the other events. 

PATROW [Ref. 3] implemented two algorithms, the time 
scale transformation algorithm and the Poisson decomposi- 
tion and gap-statistics technique, and compared them for 
computational speed and computer memory requirements. His 
results indicated that the Poisson decomposition and gap- 
statistics technique was from two to seven times faster than 
the time scale transformation algorithm, although the former 
required about thirty percent more computer memory. 

PATROW'sS work [Ref. 3] 1s also an excellent self- 
contained reference on Poisson processes, bringing many 
references together under one cover. 

A recent result of LEWIS and SHEDLER [Ref. 1] develops 
a new method for generation of points in a non-homogeneous 
Memeson process. This method, called "thinning", 1s similar 
to the general conditioning-acceptance-rejection method 
but has subtle differences which are computationally signi- 
Meeant. The thinning method is straightforward in both an 
analytical and a computational sense, and is valid for any 
type of intensity function. The thinning theorem is 
presented in Section II. 

This thesis is, in a sense, a sequel to PATROW's work 


[Ref. 3]. Its purpose is to implement the thinning algorithm 


jee 





miecomputer program form and to compare it to the Poisson 
decomposition and gap-statistics algorithm implemented by 
PATROW [Ref. 3]. The latter implementation was designed 
For a specific subset of intensity functions, degree-two 
exponential polynomials. Since the Poisson decomposition 
and gap-statistics method outperformed a general case 
megerithim (time scale transformation) by a significant 
margin, comparing the thinning method to the Poisson 
Gecomposition and gap-Sstatistics method should give a 
reasonable indication of the thinning algorithm's performance 
in generating non-homogeneous Poisson processes with other 
than degree-two exponential polynomial intensity functions. 

Section III lists the two algorithms considered, as well 
as a Special application of the thinning process which will 
be of interest to those involved in event-step simulation. 
Section IV describes the methodology used in comparing the 
algorithms while Section V deals with aspects of the thinning 
procedure which may be exploited to enhance its overall 
effectiveness in a variety of situations. Finally, Section 
VI presents the results and conclusions of the comparisons 
of the algorithms. Appendices A and B contain secondary 
Besults and computer program listings following the 


appendices. 





Pvc NNENG THEOREM 


The underlying concept of the thinning method involves 


* 
the use of a "bounding" Poisson process, ae Se > Ol 
* 
where Ny is the number of points in the bounding process in 
the interval (0,t]. This process may be either homogeneous 


ornon-homogeneous Poisson, but should be one which is easy 
to Simulate on a computer. It is called bounding because 
its intensity function, denoted ce bounds the intensity 
function A(t), of the nonhomogeneous Poisscn process which 
is to be simulated over the fixed interval (0,t']. That 


* * 
fee) = «A(t) for all t in (0,t‘]. Points at Tes 


* 


Mise, ..., Neus 


are generated for the bounding process 
over the interval (0,t']. These points are then deleted, 
or "thinned", with independent probabilities equal to 

i — (A(T, )JAa Oe Thus the Probabrlivey that a point of 


* 
the bounding process, T., iS a point of the process being 


al 
generated is equal to the ratio of the intensity functions 
Pyetuated at that point, i.e. M(T-)/A (T.). 

meme formally: 
Theorem 1. Consider the one-dimensional non-homogeneous 
Poisson process (N, eee Oe eee ees eC Meu er Lon ie. 


* 
The number of events, N Poe fhe el wecminterval  (O7t* | 


t'! 
* 
has a Poisson distribution with parameter wu (0,t') = u 


- x 
eect ) — A (0). 


1S 








* * * < 
et T., To, T 3, ee Thy" , be the times of the events 


' 
of the process in the Ree sai Ot se 

pmeeese that For 0 < t < t*', A(t) < , ae eis 
Meee 2, «6 ey Ne delete the event at T. with independent 

* * * 
Peobability 1 - A(T.) SA (T.). 

Then the remaining times form anon-homogeneous Poisson 
Mmocess With rate function A(t) in the interval (0,t']. 
BeEOOt : 

We assume that A(t) 1S continuous and use the definition 
of the Poisson process based on incremental probabilities. 
Thus we need to show that the occurrence of an event in 


(t,t+dt] is independent of the number or times of occurrence 


of events before t, and that 


PIN, way - Ny = Oo — le Hote at +eo (dt) ; 
ei Mhaseers - Nhe ~The they at. + Oldt), 

and 
SUNG ae Fe ~ Ny le ae) (ele) s 


Now we have that 


186 





mmo event from nee oem ee, erat) } 


* 
P{no event from {N, eee, Etat) i+ Pievent 


* 
from {N, we eo the (t,et+cdt) and it 1s “thinned"} 


om h(t) dt t [A (#)at]- (1 = Mey Gy. or @) (elie) 


- x # x 
nn mmc ee ett) dte— A (ty=A(t)/A (t)+dt + o (dt) 


le A(t) dt + o{dt). 


Seimilarly: 


Prone event from {N. : t OP Mii eitiet cit | 


ie 


lv 


lv 


* 
Peevent from {N, 2: t i Panis etde ewe ts not | taimmed | 


Me (e) ae -X(t)/X (t) + O(dt) 


peu) dt + oft) 


Piso it follows directly that 


Pemorerthan one event in (t,t+dt|} = o(dt) 


Moreover, an event in (t,t+dt] 1S independent of what 


happens before t because: 


dee 





* 
MeiN, : t > 0} is a Poisson process and therefore has 


independent increments, and 


2. The thinning uniform random variate is independent 


of other thinning variates, and is independent of the 


* 
Poisson process {Ny 5 5 2 Oy 
et. LD. 


megure | shows a graphical representation of a particu- 


lar case of bounding and object intensity functions. 


dle: 
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iat wGcORETHMS CONSIDERED 


foe OLooON DECOMPOSITION AND GAP-STATISTICS ALGORITHM 
Pee usage 
This algorithm is the one found by PATROW [Ref. 3] 
to be most efficient in simulation of the degree-two exponen- 
tial polynomial class of intensity functions and its imple- 
mentation by PATROW was confined to that group. Basically, 


the approach 1S to decompose the intensity function (which 


is of the form A(t) = exp (ay + a,t +f ant’) into a degree- 
One exponential polynomial function, A, (t) = exp (by + bit), 
and a difference function, Ap ft) = A(t) - A, ft). The points 


or events in the process with the degree-one exponential 
polynomial function, A, it), are generated over the interval 
Mee) Ucllizing the computationally fast gap-statistics 
erpeorithm. The points in the process with intensity func- 
eon Ap (t) are generated uSing conditioning-acceptance- 
rejection techniques. The two event streams are then 
Superposed to produce the event stream for the non-homogeneous 
Memoson process with the intensity function A(t). 

In the case where A(t) has an internal maximum 
Seeminimum in the interval (0,t'], the interval 1s par- 
titioned and treated as two separate intervals for simula- 
tion with the event streams being merged in the final step. 

Poaeronencywels Optimized by Maximizing the area 


under the degree-one exponential polynomial intensity 
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meme tion, A(t), Silegiecec=e Or COUrSEe, tO the constraint 

A, (t) Sect) 205 0 < £E =< € . This maximizes the use of 
the faster gap-statistics algorithm and minimizes the use 
of the conditioning-acceptance-rejection algorithm which 
is slow relative to the gap-Statistics algorithm. 

PATROW [Ref. 3] deals extensively with the details 
of this algorithm and it is consequently presented here 
Sey in outline form. 

2. Algorithm Statement 

meeeGactegorize the intensity function, A(t), into 

one of Six cases by examination of the coefficients a 


Jl 
and a. in A(t) = exp (ap + a,t + oer: Examples of each 
of these cases are shown in Figures 2 through Figure 7 
Mmeeated in Section IV. 
b. (1) If A(t) 18 monotone increasing or monotone 


decreasing over the interval (Cases I, II, IV and V; see 


Figures 2,3,5 and 6), decompose A(t) into Ap tt), which is 


degree-one exponential polynomial, and Ay (t) =e) — A, (t). 
Thus the decomposed functions have the forms: 
dy (t) = exp(by + bt) 
ab 2 —_ 
Ap ft) = exp (ap a,t 4p ast ) exp (bp + b,t) 
and 
A(t) = Ay ft) - Ap ft) - 


zal 





Choose Do and bs SO aS to maximize the 
area under 7, ft) Subject to Ap ft) Sp 0S)) selene sel LE ae, alg) 
et |. 


(2) Preece) “1s NOt Monotone Over the interval 


(Cases III and VI; see Figures 4 and 7), partition the 


g@nterval (0,t'] into two disjoint, contiguous subintervals, 
(0,b] and (b,t']. Choose b as the (unique) point where 
moots a Maximum (minimum) of A(t) over (0,t']. Treat 


each subinterval as in b.(1), applying subsequent steps on 
each subinterval separately, and combining results as the 
final step. 

c. Generate points in the Poisson process with 
degree-one exponential polynomial intensity function, 

A, ft), uSing the gap-statistics method. 

d. Generate and order points in the Poisson process 
Miemeintensity function Ap Ct) USinememe COnditioning— 
acceptance-rejection method. 

e. Merge (Superpose) the two event streams from 
Step 3 and Step 4. The merged stream is from the non- 


homogeneous Poisson process with intensity function A(t). 


fee tht BASIC THINNING ALGORITHM 
1. Usage 
The thinning theorem is implemented in a straight- 
£Orward manner. Typically, the bounding process used is 


* * 
homogeneous Poisson with constant rate A , where A 1s an 


Ze 





meee ound of A(t) over (0,t']. In this case efficiency 
* 

is optimized if » is the least upper bound (LUB) of X(t) 

meyer (0,t']. 


2. Algorithm Statement 


a. Generate events in the Poisson process 
(Nn, ae 7 O} with rate function heey in the fixed interval 
(0,t'J]. If the number of events generated, aa is such 
Enat aa = 0, exit; there are no events in the process 
{Ny mee > 0}. 

* * 

b. Denote the (ordered) events by Tye To re ames bas 
Set 1 = 1 and k = Q. 

c. Generate U., uniformly distributed between 0 
and l. If u; < Nae) (2 Sete mecca £0. ki) sand tT = oe 


* 
Pee ereecoual EQ itl. ir i < mn , go £0 c. 


Sk AL 


D0 One as Where n = k, and also n. 


e. Return The 


feetue ONE-AT-A-TIME THINNING ALGORITHM 
1. Usage 

In some event-step simulations, it 1S customary or 
necessary to generate only one event at a time, rather than 
an array of events. The thinning algorithm is easily 
modified to generate the next event inthe non-homogeneous 
Memescon process with intensity function A(t). In this case, 
the algorithm utilizes the time of the last event, T5 yi 
the right hand limit, t', of the fixed interval over which 
the process is being simulated; and the bounding process 


* ¢ 
intensity function, A (t). All variates are generated 
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one at a time, thus no arrays are required for storage. 
The output is Tis the time of the next event, if any, in 


the interval (T,_ eas ae 


a 


The algorithm is stated here for the case in which 


the bounding process is homogeneous Poisson, i.e., 


m (t) = a Smeelsecatmenandsan Upper bound of A(t). 
Bpeciztically: 

a is obtained by generating and cumulating 
exponential (mean = es random variates a ee wate 
Meee |, 2, ..-, UNtil for the first time, 

* x * 
Ua os S (Tey a eel pg Bi g)fA 


A detailed algorithmic statement of this procedure 
EolLilows. 


2. Algorithm Statement 


Tf i= 1, set T;_, = 0 (1-e. the left end joke) ele 
of the interval), otherwise, Ts 4 is known. Then for each 
jel, 2, ..- , the time, Tes of the event in the non- 


homogeneous Poisson process is given by the following: 


a. Set j=l 


* 
b. Generate E. 4 an exponential random variate 
P 
k 5S) 
feeea mean 1/A . If T. + E. is greater than t', 
hb eet nly 
exit; there are no more points in the interval (T,_,/t'l- 
c. Generate Clay uniformly distributed between 0 
mac lL. if 
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set 


and exit. 
See Ochnerwise set j = jt+tl; go to b. 
Note: U; ; and o are uniformly distributed 
t 


between O and l. 
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IV. METHODOLOGY FOR ALGORITHM COMPARISON 


ew eAoURES OF EFFECTIVENESS 

Two quantifiable measures of effectiveness were chosen 
as yardsticks for algorithm comparison. These were compu- 
tational speed and computer memory requirements. Some other 
considerations, such as programming ease and robustness, 
are discussed in Section VI. It must also be recalled that 
the classes of intensity functions for which the two 
algorithms are usable are different. The Poisson decomposi- 
tion and gap-statistics algorithm is only easily implemented 
for a restricted set of intensity functions, those of the 
form A(t) = exp (ay + a,t + sede 1.e. the degree-two 
exponential polynomial. Conversely, the thinning algorithm 
1s valid for any positive intensity function. Thus a direct 
comparison can only be made in that subset of intensity 
functions for which both algorithm implementations are 
valid, the degree-two exponential polynomials. 

PATROW [Ref. 3] developed six sample intensity functions, 
all special cases of the degree-two exponential polynomial, 
and these are used herein as the test cases. These are 
described in Section IV.C.3 below. 

1. Computational Speed 

Typically, computer time is a costly commodity 
mieeconomic terms. It may also be a significant factor 


in determining more mundane considerations such as job 
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Priority and thus turn-around time. Thus computer run 
time is a natural candidate as a measure of effectiveness 
for comparing competing algorithms. 

PATROW 9[ReE. Si) utllazed a procedure in which 
event streams from each of six sample intensity functions 
were replicated several times in "packages". The number 
of replications was large if the expected number of events 
in the event stream was small, and vice versa. Thus the 
product of the number of events times the number of replica- 
tions was kept on the same order of magnitude. For simplicity, 
the same technique was used here although results showed a 
wide variation in the run times for the six packages. 

Programming of the thinning algorithm was done so as 
to minimize run time while maintaining parity with the 
Poisson decomposition and gap-statistics algorithm wherever 
direct comparison could be made. For example, shuffled 
random numbers are called in both programs and both are 
dimensioned to accommodate event streams of up to 5000 events. 

Undoubtedly further programming refinements exist 
which might increase slightly the speed of one or the other 
algorithm. Also, different computers might have unique 
features which could be exploited. The overall purpose here 
was to obtain a relative order of magnitude comparison and 
it believed that this objective was accomplished in every 


meaningful sense. 
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2. Computer Memory Requirements 

This is the second obvious means of comparing two 
algorithms. Again, some core reduction could undoubtedly 
be made by a sSsophiSticated programmer. Most notably, core 
requirements can be reduced substantially if only one non- 
homogeneous Poisson variate is generated at a time (the 
one-at-a-time algorithm) but this has the predictable 
effect of increasing execution time considerably (see 


Tables IV, V, and VI). 


B> MEASUREMENT CONSIDERATIONS 

Measurement of computer memory requirements is straight- 
forward and deterministic. 

Measurement of computational speed, more specifically 
Central ProcesSing Unit (CPU) time, is quite another matter. 
Pest Of all, the number of events in each replication of 
the nom-homogeneous Poisson process varies causing CPU time 
to be a random variable. More important, however, are the 
effects of internal computer procedures. 

In the first place, the so-called CPU time printed out 
on the normal IBM-360 output has only a general relationship 
to the actual computational time required by the CPU. This 
is caused by the addition of certain "overhead" time. This 
overhead time is a function of the number of other programs 
in the system as well as such factors as compilation and 
linkage times. Thus the same program run at two different 


times may differ in "CPU time” by a factor of two or more. 





Program execution times were isolated from compilation 
and linkage times by the use of a system subroutine, 

GETIME. This subroutine allows the user to initialize an 
internal timer within the program and read cumulative time 
at various points in the program. Although this method is 
not exact, it does measure actual elapsed CPU time to within 
a small fraction of a second. This does not, however, 
entirely alleviate the time-of-day effect experienced when 
running the same program at different times. That is, 
although the elapsed CPU time can be measured accurately, 
the same program will generally have somewhat different 
execution times each time it is run [Ref. 4]. Theoretically, 
the execution times would be constant for stand-alone runs, 
1.@€. runs with no other competing programs in the system. 
This is rarely realized in practice. 

These considerations lead to the development and use of 
the side-by-side setup described below. This method appears 
to be statistically sound as a means of dealing with the 
problems of time measurement. Due to the differences in 
execution times noted, the best measure of effectiveness 
was determined to be a ratio of execution times for the 


respective algorithms, rather than absolute times. 


See LBoL SETUP 
1. Computational Speed 
The central idea here was to equalize the effects 
of non-essential processes on each algorithm. This was 


accomplished by the following algorithm: 
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ime cet k = l. 

2. Zero internal time clock. 

Seeeecall Algorithm A. Replicate M times. 
4. Read internal time clock. Store time. 
5. Zero internal time clock. 


Seeeecall Algorithm B. Replicate M times. 


7. Read internal time clock. Store time. 
Pemeoec kK = k+il. If k is greater than aert Gone o. 
Otherwise go to 2. 
9. Compute mean and variance of the Lane execution times 
fmemeeeach algorithm. 
10. Compute ratio of means. 
This procedure was used in all comparisons. M, 


the number of replications per package, varied between 30 
and 100 as discussed above. res eye the number of times each 


package was replicated, was typically set equal to thirty. 

2. Computer Memory Requirements 

To measure computer memory requirements, a small 

Main program was written, calling the subroutine which 
implemented the program being measured. Total program length 
in bytes was obtained from the standard computer output and 
the core alloted to the main program was subtracted to 
obtain the desired figure. This includes all library rou- 
tines and arrays for storage of event times and arrays of 


Fandom variates. 


30 





The core requirements are deterministic in that 
they do not change from one run to another but are strictly 
Peerumetion of the program coding. 

3. Test Cases Utilized 

PATROW [Ref. 3] developed six sample intensity 
functions representing the possible variations in sign and 
relative magnitude of the coefficients ay and a, in the 
Meeonential polynomial, exp(a, + a,t + age) 

Since the sample intensity functions were designed 
to test different aspects of the Poisson decomposition 
and gap-statistics algorithm, they were also used for com- 
parison here. Although each algorithm is affected by 
different considerations, the test cases do, coincidentally, 
Me@eethne thinning algorithm through its paces. 

For continuity, the six test cases, or sample 


intensity functions, are presented below in Figures 2 through 


ge 


om 
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V. EFFICIENCY AND PROGRAMMING CONSIDERATIONS 


A. GENERAL 

This section deals with factors which affect the per- 
formance of the thinning algorithm. Four specific areas 
are presented in which significant gains in terms of com- 
putational speed may be realized. With the exception of 
section V.D (which applies only to the case of exponential 
polynomial intensity functions) these considerations apply 
to the general class of intensity functions. 

maegeneral application, one of the primary indicators 
of efficiency is the relative size of the area under the 
Mieenisity function to that of the bounding function, i.e. 


jm@enratlio, R, given by: 


fee og 
* 
Ree ee (Ede fa Note) dt 
0 0 


Since both numerator and denominator are simply the respec- 
tive integrated rate functions, A(t), evaluated at either 
end of the interval, Ris the ratio of the expected number 
of events in the two processes, l.e. E(N,,1/E(N,,1- 

Case 1 of the sample intensity functions is particularly 
gilustrative (see Figure 2). The intensity function 


mee = exp(1.6 + 0.015t + 0.0005t*) 1s bounded on the 
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interval (0,100] by a least upper bound (LUB) of 3294.47. 
If a homogeneous Poisson process with rate equal to the 
Meets used as the bounding function, BIN, 1] = 329,447 
points will be generated on the average. Of these, all 


but 1464 will be rejected on the average (i.e. E[N,,]). 


t! 
The ratio of the respective expected values is thus 
m7 329,447 = 0.0044 = the ratio of the areas under the 
intensity functions. 

Thus a rough relative measure of the efficiency of the 
thinning process in a particular situation can be gained by 
examining a graph of the two intensity functions, even if 
the expected values are not easily calculated. This pro- 
cedure may also be an indicator in deciding whether to 


partition the interval and use different bounding functions 


Smmeach subinterval. 


B. UTILIZATION OF ARRAYS OF RANDOM VARIATES 


Computer generated random variates are used both in 
* 


Semerating the points of the bounding process ve 


On 
m@a@ein the actual thinning process itself. Since the 
number of variates required is typically large, efficient 
generation becomes a programming consideration for medium 
to large scale simulations. 

The basic thinning program presented herein requires 
both exponential and uniform random variates. Both types 


are obtained utilizing the random number package, LLRANDOM, 


developed by LEARMONTH and LEWIS [Ref. 5]. Shuffled 


BE. 





exponential random variates of mean 1.0 are generated uSing 
the SEXPON subroutine while shuffled uniform (0,1) random 
variates are obtained from the SRAND subroutine. Both of 
these routines offer considerable "economies of scale" in 
terms of time when multiple numbers or arrays of variates 
are generated at once, aS opposed to one-at-a-time genera- 
tion. Using the test setup of Section IV.C, average times 
to generate varying quantities of random numbers were 
determined. Table I reveals the relative savings realized 
by calling large arrays of random numbers. Thus considerable 
time can be saved by generating ail required random numbers 
meereone Subroutine call. Programming difficulties involve 
deciding how many variates to generate. The general goal 
is to generate as many aS needed while keeping the unused 
excess to a minimum. The balance used was to generate the 
expected number required plus an excess of four standard 
deviations. For example, in the generation of the bounding 
process, the expected number of points, EIN, 1] 1s eee! 
and the variance is the same. Thus the number of exponen- 
tials called was y + 4Vy where y = \ oedy Provision is 
Made for the unlikely (1 in 40,000) case that more are 
required. 

For specific applications this procedure could be 
improved slightly. For example, if the expected number of 
Bounts, EIN, J, is small, e.g. 100, then the expected 


excess (four standard deviations) comprises forty percent 
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Type of Number Total Time Mean Time Per 


Variate Called (usec) Variate (usec) 
Exponential i 784 784 
Exponential 10 1293 WA 
Exponential 100 7343 72 
Exponential 1000 68046 68 
Una torm 1 2 1.3 2a: 
Pama torm 10 eis 1. Lae 
Uniform 100 B2n6 Se 
Uniform 1000 21544 Ze 


Sample Size = 200 (each grouping) 


Table I 


Generation Times for Arrays of Shuffled Random Variates 
From LLRANDOM 


* 
of the total whereas for large E(N,, 


expected "waste" is only about six percent. In the former 


], e.g. 4000, the 


case, reducing the "padding" to one or two standard devia- 
tions would, on the average, increase efficiency slightly 
Peeenough the probability of a second subroutine call for 
more random variates would increase. 

AS an example, if we were to call 100 exponential 
Variates one at a time, the total time is 78,400 usec, 
compared to 7,343 usec for 100 exponential variates called 
im an array. For 100 + 4Y100 = 140 variates, the time 


is still small compared to 78,400 usec. 
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feo ri LIZATION OF INTENSITY FUNCTION LOWER BOUND 

One of the most time conSuming repetitive operations 
is the computation of the intensity function value, A(t), 
during the thinning process. In the case of the exponential 
polynomial intensity function, this involves one power, two 
multiplications, two additions, and one exponentiation for 
each point generated in the bounding process. Since points 
are accepted forthe nomhomogeneous Poisson process when- 
ever the uniform (0,1) thinning variate is less than the 
ea 1.0 MED /) (t), considerable time Savings result if the 
Mieensity function has a positive lower bound, Say A, 
Since points are always accepted when the uniform (0,1) 
Variate is less than the ratio View In the general case, 
this ratio must be calculated only once. The expected num- 
ber of intensity function computations which are alleviated 
by the use of the lower bound is given by (AJ) EIN, where 
A is a lower bound of the intensity function; “ 1S an upper 
bound of the intensity function (both bounds are over the 
interval (0,t‘']) and EIN, 1] is the average number of points 
to be thinned, i.e. the average number of events in the 
bounding process. 

eos Clear that the cloSer 4 is to being the greatest 
lower bound (GLB) and the closer * is to being the LUB, 
the more efficient the program. 

If the intensity function is strictly non-decreasing 
a further (and potentially great) improvement is realized 


by initially setting A equal to A, and then setting it 
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subsequently equal to the last value of the intensity 
function, A(t). This results in a monotone increasing 
lower bound and thus a decreasing probability of evaluating 
the intensity function. 

Test cases II through VI were run Side by side with and 
without the use of a lower bound for the intensity function. 
On the average, the program which did not utilize a lower 
bound required twenty percent more time than the program 
using a lower bound. Please see Appendix B for case-by- 
case comparison. 

D. UTILIZATION OF EXPONENTIAL VARIATES FOR THINNING OF 

EXPONENTIAL POLYNOMIAL INTENSITY FUNCTIONS 

The time requirements for evaluating A(t) were discussed 
in Section V.C above. In the case of exponential polynomial 
Mimeensity functions, e.g. A(t) = exp (a, at a,t + age): the 
major contributor to computation time is the exponentiation 
Operation. Exponentiation can be avoided by utilizing the 


following relationship: 


* 
U, < heey. ibiestgiel Chale alae 


* * ? 
-In U; > ink - mwece) (= ini ss (a, t+ a,jt + aot), 


gs) 
ll 


where: 


U. icomcwuineeeonrm (0,2) randemevariate, 


E. is an exponential random variate with mean one. 
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Maus the thinning test to accept points from the 


bounding process becomes: 


* 2 * 
‘ = a + = ; 
jiag E. > Ina (ay a,t + aot ymogece T5, accept 
* 
T, as a point inthe non-homogeneous Poisson process 


* 
with rate A(t); otherwise, reject T; (pi eeret ce le. Ante): 


The key to this relationship lies in the fact that if 
Mmemamstributed uniform (0,1), then -ln U is distributed 
as a unit exponential variate, i.e. an exponential variate 


with mean one. This is shown by the following: 


meee be uniform (0,1). 


MeperiU < x} = P{lin U< In x} = P{-ln U > -ln x} 
Peer iu < x} = x, thus let y = -ln x, 
Maen Pi-ln U > y} = exp(-y). 


Thus -ln U is distributed as a unit exponential variate. 


Although more time is required to generate the exponential 
memagom variates for thinning than the uniforms, the alle- 
Miaeron of the exponentiation operation more than compen- 
Sates for the additional generation time. This is because 
SEXPON, the portion of LLRANDOM which generates exponentials, 
generates exponential variates by the Marsaglia "rectangle- 
wedge-triangle" method, which is faster than taking logarithms. 

Since exponential random variates are used in the genera- 
tion of the bounding homogeneous Poisson process, an addi- 


tional time savings can be realized by using the variates 
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Mmmem are left over (i.e. not used) from generating the 
bounding process (these are generated in arrays). 

For the test cases considered, use of exponentials 
for thinning resulted in an average time savings of ten 


percent. Please see Appendix B for case-by-case results. 


feecacYyCLING OF THINNING VARIATES 

As mentioned above, a uniform or exponential random 
variate iS required for each point to be "thinned". Each 
of these variates requires a Significant amount of time for 
generation. Obviously a time savings would be realized if 
fewer variates were required. 

me Recycling of Uniform Variates 

Assume U. mmr rorm Or) ebue that LES value is 

unknown. Assume then that further information becomes 
available that U; eS Sac oe) eSue tts Value 
is still unknown. Then U. is uniformly distributed over 
Mme interval (0,a). If U. 1s then computed by “scaling 


1i+l1 


sol Us, ioe .aividing U. by a, then Us, WS Wee Orme Oka) . 


Hl 
Semiiarly, if U. is uniform (0,1) and subsequent information 


places it somewhere above a, then Us, ~ (U, - a)/(l - a) 


il 
Memmi torm (0,1). Thus by conditioning on whether the 
variate is greater than or less than a given value, a new 
variate can be computed with the desired properties. 
Moreover, this variate is independent of its predecessor. 


Pe oneowtntinmang algorithm, each point 1s tested 


* ae 
meg a uniform (9,1) variate. Specifically, if U. S A(T.) /A 
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the point i is accepted as a point in the non-homogeneous 
Poisson process. Since the ratio A(T.) /A- is between zero 
and one, and the only test is whether U; 1s less than or 
greater than the ratio MEIN the next uniform (0,1) 
variate, Us ay! can be generated uSing the rules above. 


The algorithm is: 


ioe 6CLet U; Demuna corm (Olja Lt U, is less than 
* * * * 
a= A(T, )/ Palet Uiay = U,/ (A(T) / i exit. 
* * x Ok 
2. Otherwise let Us4y = LU, ~ (ACT) 7% IG at ele) Val 
Us4y Menumicorm (0,1). 
In theory, only one uniform random variate is 
Meeered for the entire thinning process. In computational 


practice, however, care must be exercised because of the 
finite capacity of the computer to represent numbers. After 
ten to twenty divisions the scaled uniform number will 
consist only of low-order bits of the random number and 
these are usually not uniformly distributed. 

If the intensity function has a poSitive lower bound, 
further efficiencies can be gained, in combination with the 
Meoceacures of Section V.C above. Since multiplication is 
computationally faster than division, the value L/(A/A) = UN 
can be precomputed and stored. Thus if U, < Wile 
Gray ~ cen can be computed as the next thinning variate. 
Meee that no intensity function calculation is required. 
However, if U. > Ve the thinning method proceeds with 
the next step, evaluating the intensity function at the 


* « 
point T. and determining whether or not to thin the point. 


46 





Now, further information is known about U.. Specliheal ly, 
; * * ; x ; 
either U. > A(T.) , in which case T. LS EnanmMed Or 


* * 


* 
ee J; < AtT,)/A . In either case le 


é can be computed 


1 
Pye scaling up” U.. 

iio, tme algorithm for recycling uniform random 
Mertates for thinning is as follows: 


* * 
Beet U. < A/d , let U,,, = U,-A /A and exit. 


* * = * x 
Semmes A/A 6< Us < ACT.) fr Pole tr Us44 = (i "U,-A)S (A(T; ) -A) 
and exit. 
* * 
3. Otherwise, U. > MT.) /A PaLet 


a x * x x 
U meee J, WACT.))/(A ACT.) ). 


i+] 

By precomputing MOM, PoromEeecve ling Drocedure 
requires only one multiplication in the case where 
U; Ss Dy . Otherwise one multiplication, two subtractions 
and one division are required. In either case the recycling 
procedure is generally faster than generating uniform 
random variates from a random number generator, even when 
a logical IF statement is added to check for extreme 
values ("small bits"). 

mmeesccycling Of Exponential Variates 
This section applies only where the intensity func- 

tion is exponential polynomial. Here the possibilities 


mee less promising. In the general case where no lower 


meund, 4 is used, the following algorithm would apply: 


* ~ 3 * 
if E; Ean Ae in ACT), let Bead = E - In X} + 1n A(T, ) 
Otherwise 

_* * * * 
Let E,,, = In (AX - A(T, VSO ‘exp (-E,) - A(T, )) 


where E. is a unit exponential random variate. 


A7 





* CaF 
In the first case, E. Beem A =) in A(T), a time 
x 
savings would generally be realized since In X could be 


* 
computed once and stored and ln A(T, ) 1S Simply the value 


* * 
of the polynomial, i.e. Ao aE a,T; + aa which must be 
computed in any event. In the second case, however, the 
cure is truly worse than the illness. It is faster to 


Simply generate another exponential variate, assuming they 
are called in arrays. 

For there to be a time savings, however, it must 
be possible to make a reasonable prediction of the number 
of exponentials which must be generated. Otherwise an 
excesSive number of calls to the random number generation 
Subroutine may destroy the gains made through recycling. 

Gimerc case where a lower bound, A, for the inten- 
eevee tunction exists and is positive, it 1S possible to 


determine the expected number of exponentials which must 


be generated. Variates are reused if they are greater than 

* * * 
Mme 7). That is, if E, > ini), then Bea. > E, - NENA 7 A) 
Otherwise, a new (1.e. non-recycled) variate is used. Thus 


mies orobability of not recycling, p, is wee and the number 
of variates required is binomial with mean np, where a 
1s the number of points to be thinned. 

Empirical results for the five test cases considered 
are shown in Appendix B. Using the calling rule of expected 


number plus four standard deviations for generating thinning 


exponentials yielded inconclusive results as compared with 
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the procedure in which exponential thinning variates are 
generated in arrays with no recycling. As expected, for 


larger N (Cases Il, Tif and V), recycling provided a slight 


i 
time advantage (seven percentage maximum) while for small 


N (Cases IV and VI) recycling was slower. In case VI, 


t! 
recycling caused run time to be approximately five percen- 
meee greater than that without recycling. Using a calling 
rule of expected number plus two standard deviations reduced 
the disadvantage slightly to four percent. The reason that 
recycling can cause longer run times than not recycling is 
that an additional logical IF statement is required for the 
recycling program. Again, when exponentials are used for 
thinning and the mean number of points to be thinned is 
Smmeene Order of two or three hundred, it is probably not 
worth the effort to recycle. If Several thousand "thinnings" 
are required, the savings may indeed be worthwhile. 

Results were somewhat surpriSing in the general 
class of intensity functions for which uniform variates are 
used for thinning. It was expected that a significant savings 
would be realized since the uniform variates can be recycled 
memati Situations. In fact, test runs were run in which 
only one uniform variate was generated for the entire run 
fee recycling used throughout. The only program statement 
which added time was a logical IF statement to preclude 
mavading by zero and to diminish the probability of small 
bit usage. As expected, some bias was experienced in the 


mean and an unusually high variance was noted, indicating a 
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low degree of "fidelity" to the true non-homogeneous Poisson 
process being simulated. Of particular interest however, 
were the results on execution time. Since this setup 
essentially gives the lower bound on execution time for 
recycling of uniform variates. it was expected that signi- 
ficant time savings would be realized in comparison to no 
recycling. In practice, however, the savings were minimal, 
with a maximum of only three percent savings. This is 
attributed to the efficiency of the LLRANDOM package in 
generating uniform variates with the logical IF statement 
being only a secondary cause to longer execution time. 

Appendix B shows results for both exponential and 
moeeorm thinning variates. 

The key point to keep in mind here is that the above 
results reflect the case where thinning variates are called 
in arrays, i.e. many at a time. Thus the comparison is 
between a very "fast" variate and recycling. As discussed 
above, calling by arrays results in considerable time savings 
compared to one-at-a-time generation (up to fifty times 
faster). Thus, in the case where a slower random number 
generator is utilized or where variates are called one at 
a time, use of the lower bound may indeed result in con- 


Siderable time savings. 


FF. FINAL PROGRAM 
A general program was developed which incorporated the 


efficiency considerations discussed above. The program is 
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general in that it can be used with any of the general 
class of intensity functions, whether exponential polynomial 
Or not. The program is essentially four programs, each 
used in a specific case. The program classifies simulations 
into the four classes by asking two questions: 

1. Is the intensity function exponential polynomial? 

2. Does the intensity function have a positive lower bound? 

The first of these determines whether uniform variates 
(general case) or exponential variates (exponential poly- 
nomial case) are used for thinning. The second consideration 
merely deletes an unnecessary logical IF statement in the 
case where no lower bound is used. 

The computer program, NHPP, is listed after the appen- 
dices, and requires a user supplied subprogram FUNCTION FCN(T) 
to compute the intensity function values for each value of 
feeeeele the intensity function is exponential polynomial, 
only the exponent portion should be calculated, i.e. the 
statement FCN = ay + a,t Ze ay (for degree-two polynomial) 
should appear in the subprogram. Otherwise the entire 
intensity function should be evaluated. 

Empirical results showed that the final program, 
utilizing the efficiency considerations mentioned in this 
Mmepter, resulted in a program which ran in two-thirds the 
time of the basic thinning program. 


Please see Appendix B for case-by-case results. 
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VI. RESULTS, CONCLUSIONS AND RECOMMENDATIONS 


A. GENERAL 

This section presents the results of comparison of 
the Poisson decomposition and gap-statistics algorithm 
with three variations of the thinning algorithm. These 
Variations are the basic thinning algorithm, the modified 
thinning algorithm (final program) and the special case 
one-at-a-time algorithm. Section B presents the performance 
of each of the algorithms when measured by the two measures 
of effectiveness, computational speed and computer memory 
requirements. Section C examines the results with a view 
toward identifying the strong and weak points of each 
algorithm. Section D recommends further avenues of study. 

Again, in comparing the two classes of algorithms, one 
ere distinction must be kept in mind. That is that the 
Poisson decomposition and gap-statistics algorithm as imple- 
mented by PATROW [Ref. 3] is limited to a special class of 
Micensity functions, i.e. exponential polynomial of degree 
mye (or less). Although the algorithm could be adapted to 
Meegmer Order polynomials (by further bisection of inter- 
vals), the already complex programming considerations would 
meow Significantly. In contrast, the thinning algorithm 
is a completely general method which is analytically valid 
for any functional form of permissible intensity functions 


Meeottive and right continuous). 


or 





The results presented here are necessarily limited to 
mieemelassS Of intensity functions for which both algorithms 
can be compared, i.e. the degree-two exponential polynomial 
class. The purpose was to determine the relative performance 
Of the thinning algorithm on this piece of common turf with 
the heretofore champion, Poisson decomposition and gap- 
meacistics. 

Mr basic result is that the thinning algorithm is 
indeed quite competitive with the Poisson decomposition and 
Sepeseacistics algorithm in the area of mutual validity. 
This, combined with its ease of programming and ability to 
generate variates from any intensity function, make the 
thinning algorithm a highly attractive tool for generating 
non-homogeneous Poissonprocesses of any type. 

One shortcoming of the thinning program was revealed by 
the first test case considered (see Section IV.C). This is 
a fast rising exponential polynomial which rises from a 
Value of five to almost 3300 over the interval (0,100]. For 
Mi = 3294.47, the expected number of points in the bounding 
process is 329,447 while the non-homogeneous Poisson process 
being simulated has an expected number of only 1464 points. 
Thus all but one point in about 200 are thinned out. The 
thinning algorithm could be more efficiently adapted to 
this case by partitioning the interval, alleviating the 
necessity to store over 300,000 bounding process points. 


However, the efficiency involved would still be low, and the 
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best solution appears to be to utilize the Poisson 
decomposition and gap-statistics approach. The key point 
here is that the problem is easily recognized beforehand, 
as discussed in Section V.A, and avoidable. 

Table II presents a general comparison of the two 


algorithms. 


Seo URES OF EFFECTIVENESS RESULTS 

Chapter four details the comparison procedure utilized 
to develop the following results. 

ime ine Basic Thinning Algorithm 

Salient features for this case include the use of 
uniform variates for thinning and one-at-a-time generation 
of exponential and uniform variates. 

Table III presents the results for each of the five 
test cases run. Algorithm A is computer program DEGTWO, 
miemr Olsson decomposition and gap-statistics program 
developed by PATROW and listed at the end of this paper. 
Algorithm B is computer program NHPTHN, the basic thinning 
algorithm, also listed. 

The thinning algorithm was fastest in two out of 
the five test cases run and required eighty percent of 
the core space required for the gap-Statistics algorithm. 

Table VII lists core storage requirements for each 
aeroorithm. 

fee tne Modified Thinning Algorithm (Final Program) 
This section compares the best case performance 


mer both algorithms. 
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The modified thinning algorithm includes: 

1. Use of exponentials for thinning of exponential 
polynomial intensity functions 

2. Use of lower bounds 

3. Partial recycling of exponential thinning variates 

4. Use of exponential variates left over from generation 
of the bounding process. 
Fach of these refinements are discussed in Section V. 

The Poisson decomposition and gap-statistics 
algorithm used was again the implementation by PATROW 
Meets 3], program DEGTWO. In addition to the normal running 
of the program, a second set of comparisons was made uti- 
lizing separately calculated values for ae Ene Donic. FOr 
the conditioning-acceptance-rejection routine. This is 
discussed by PATROW [Ref. 3]. These runs are indicated 
by aSterisk (*) in Table IV. 

In the first case, the thinning algorithm was faster 
in four out of the five test cases, with the best relative 
Be-erormance occurring in Case VI. 

When the improved values of = were incorporated 
Poeecases IV, V, and VI, the Poisson decomposition and gap- 
Statistics algorithm improved substantially, winning in 
three out of the five test cases. 

The relative advantage in computational speed was 
less than a factor of two in all cases with a maximum of 


ms5ecol in tcavor of the thinning algorithm in case VI. 
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Core storage for the thinning algorithm was determined 
by using only the part of the algorithm which used exponen- 
tial thinning variates. This precluded the requirement for 
storing 5000 uniform variates which are not used. The 
Poisson decomposition and gap-statistics program (DEGTWO) 
required about 88,000 bytes of core storage as compared to 
feeomes+,000 for the thinning program (NHPP modified to 
exclude unused uniform variate array). Detailed results 
are listed in Tables IV and VI. 

fee the One-At-A-Time Thinning Algorithm 

As discussed in Section II.C, the one-at-a-time 
algorithm was developed only to test the relative efficiency 
of the algorithm used to generate the next event in a non- 
homogeneous Poisson process. This latter requirement may 
arise in event-step Simulation where only the next event 
in anon-homogeneous Poisson process is desired rather than 
an array containing all events in a specified interval. 

Computationally, the one-at-a-time algorithm is 
ere Similar to the basic thinning algorithm. The only 
essential difference is that the basic thinning algorithm 
generates and stores all the points in the bounding process 
femcensity function ee before thinning, whereas the one- 
at-a-time algorithm generates a point in the bounding process 
and thins that point before continuing. The latter method 
removes the requirement for an additional array to store 


the bounding process points. This in turn saves about 


2S), 





20,000 bytes of core storage requirement when the programs 
are dimensioned to accept 5000 points. 

As implemented here, the one-at-a-time algorithm 
Simply generates the next point in the non-homogeneous Poisson 
process and stores it, stopping when the last point generated 
lies outside the interval. All variates in this program are 
generated one at a time. The results shown in Table V are 
thus a good indicator of the relative efficiency of using 
this method. As can be seen, the one-at-a-time algorithm 
(program NHPOAT) is faster than the Poisson decomposition 
and gap-statistics algorithm (program DEGTWO) in three of 
the five test cases run. This 1S true despite the fact that 
DEGTWO generates all variates in arrays, taking advantage 
of the time economies of scale mentioned in Section V. The 
One-at-a-time algorithm also requires forty percent less 
core space. 

From the tables one can also see that since both the 
best case thinning algorithm (Table IV) and the one-at-a- 
time thinning algorithm (Table V) are compared to the 
Pemssonm Aecomposition and gap-statistics algorithm, it is 
possible to obtain a reasonable comparison of the best case 
thinning algorithm and the one-at-a-time thinning algorithm. 
For example, for the sample intensity function used in 
Mase ff (see Pigure 3), the ratio (.8127/.5541) = 1.47 
indicates that execution time for the one-at-a-time thinning 


algorithm is almost fifty percent greater than that of the 
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Beoe case thinning algorithm. The ratios for Cases III, 
emmeestia Vil are 262, 1.27, 1.37, and 1.21 respectively. 

The tables also demonstrate the time-of-day effect 
discussed in Section IV. That is, the execution time for 
a given program is not the same each time it is run. For 
example, program DEGTWO took 21.8808 seconds for the run 
recorded in Table V compared to 18.8917 seconds for the run 
recorded in Table IV. For this reason, ratios of execution 
times were chosen as the measure of effectiveness rather 
than absolute times. 

A FORTRAN program listing, NHPNXT, is provided at 
the end of this thesis. This program generates the next 
point ina non-homogeneous Poisson process with a user 
Supplied intensity subprogram, FUNCTION FCN. This program 
can be uSed in conjunction with event-step simulation 
programs, including SIMSCRIPT, where it is desired to mini- 
mize core space (at the expense of speed) or where only one 
event is desired at a time. Core requirements are shown 


maetaple VI. 


eee CONCLUSIONS 

Both the thinning algorithm and the Poisson decomposition 
and gap-statistics algorithm include two general types of 
Operations: a "generating" process and a "Second stage". 
For the Poisson decomposition and gap-statistics algorithm, 
the generating process is the non-homogeneous Poisson process 


with degree-one exponential polynomial intensity function. 
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For the thinning algorithm (as implemented herein) the 
generating process is homogeneous Poisson. 

The second stage for the Poisson decomposition and gap- 
Statistics algorithm is the actual decomposition and genera- 
tion of variates by the conditioning-acceptance-rejection 
method. For the thinning algorithm, the second stage con- 
Sists of the thinning of the points in the bounding process. 
Thus one algorithm generates events and adds more events 
from a second process while the other generates events and 
Subtracts some out. 

The strongest point in the Poisson decomposition and 
gap-statistics algorithm is the highly efficient generation 
of the events in the degree-one exponential polynomial inten- 
Peeve cunction process. This is done with the gap-statistics 
peigorrcthm which is two to five times faster than the thinning 
algorithm for this type of process (see Appendix A). At 
the same time, the conditioning-acceptance-rejection routine 
is relatively quite inefficient. 

There are many considerations in predicting the relative 
Success of the two algorithms, i.e. which will be faster in 
a given situation. For example, the Poisson decomposition 
and gap-statistics algorithm is affected by factors such 
as whether or not partitioning is required, the percentage 
of the total number of variates which come from the degree- 
One exponential polynomial process, and whether time rever- 


sal is required. For the thinning algorithm, the fraction 
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of the lower bound divided by the upper bound Wa would 
seem to be a good indicator of success. 

For the test cases considered, however, the only con- 
Sistent indicator was the expected number of events in the non- 
memogeneous Poisson process being simulated. The smaller the number 
of events inthe non-homogeneous Poisson process being simu- 
lated, the better the relative performance of the thinning 
algorithm over the Poisson decomposition and gap-statistics 
algorithm. Thus it appears that each algorithm has a fixed 
and variable part in terms of time. The thinning algorithm 
has a shorter "setup" cost in terms of time but the variable 


cost or "cost per additional variate" seems to be smaller 
for the Poisson decomposition and gap-statistics algorithm. 
The exact cause of this phenomenon is not known although 
it appears to be centered in the conditioning-acceptance- 
rejection routine. 
In the larger spectrum of non-homogeneous Poisson process 


generation, it seems clear that the thinning algorithm is 


the best all-around method available. 


D. RECOMMENDATIONS 
Two specific areas for further study are recommended. 
First, the thinning algorithm as implemented here uses 
only homogeneous Poisson processes for bounding. It might 
be worthwhile to investigate the possibility of using other 
processes, such asnon- homogeneous Poisson processes with 


degree-one exponential polynomial intensity functions, as 
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bounding processes. This would allow the efficient gap- 
Seaeastics algorithm to be utilized although function 
evaluation would in general become more time consuming. 

The second area is in finding the optimum method for 
generating the degree-two exponential polynomial class of 
intensity function. These will undoubtedly remain of 
interest due to their statistical properties. Here, it 
seems clear that the best features of the two algorithms 
can be combined. Specifically, the Poisson decomposition 
and gap-statistics algorithm can be modified to use thinning 
rather than conditioning-acceptance-rejection for generating 
the points in the difference function process, i.e. the 
process with intensity function ds =e) = A, (t). Also 
fiemweriterion for the decomposition might preferably be 
that the intensity function of the remainder be monotonically 
increasing. This would make it easy to find the upper 
bound for the function, and the most efficient version of 
the thinning algorithm, that where the number of computations 


of the intensity function is minimized, could be used. 


66 





APPENDIX A 
GENERATION OF DEGREE-ONE EXPONENTIAL 
POLYNOMIAL INTENSITY FUNCTIONS 

The generating process for the Poisson decomposition 
and gap-statistics algorithm is anon-homogeneous Poisson 
process with degree-one exponential polynomial intensity 
function. This is generated by using the gap-statistics 
method, which is subroutine NHPP2 in the DEGTWO program 
(see listing below). 

To determine the relative speed of the thinning algorithm 
compared to the gap-statistics algorithm, two simple degree- 
one exponential polynomial intensity functions were developed 


and simulated. 


Table VII presents the results. Case A is a monotone 
decreasing intensity function, A(t) = exp(3.4 - 0.02°-t) 
over the interval (0,100]. Case B iS a monotone increasing 
Mm@eemsity function, A(t) = exp(.693 + 0.03-t) over the 


Ga@eecvyal (0,50). 
Results show that the gap-statistics algorithm is from 
two and a half (Case B) to four and a half (Case A) times 


faster than the thinning algorithm. 
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APPENDIX B 


RESO On meer NCy MODIFICATIONS 


This appendix presents, in tabular form, the results 
of comparison of the programming modifications listed in 
Section V. 

Table VIII shows the effects of utilizing lower bounds 


meteeene intensity function. Test conditions include: 


IT. Use of uniform thinning variates 
2. Recycling of thinning variates 
3. Use of arrays of variates 


Table IX shows the gains realized by employing exponen- 
Mmeebecninning variates in contrast to uniform thinning 
variates for exponential polynomial intensity functions. 


Teoe conditions include: 


1. Use of arrays of variates 
fee Recycling of thinning variates 
Peeeeuse Of lower bounds for intensity function 


Table X shows the results of recycling versus no 
recycling where uniform variates are used for thinning 
while Table XI shows the same comparison when exponential 
variates are used for thinning. For both cases, test 
conditions include: 

1. Use of arrays of random variates 


2. Use of lower bounds for intensity function 
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Table XII presents the results of incorporating all 


of the programming improvements into program NHPP. The 


final thinning program, NHPP, 


1S compared to the basic 


thinning program without modifications, NHPTHN. The 


Pesential differences are: 


NHPP (final program) 


Arrays of variates generated 
Exponential variates used 
maemechinning 


Lower bound of intensity 
met 1OnS utilized 


Thinning variates recycled 
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NHPTHN (basic program) 


Variates generated one 
at a time 


Uniform varlates used 
cone | glqnemmate, 


Lower bound = 0.0 


No recycling used 
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